#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include "galois_add.h"
#include "galois.h"

static unsigned char *table_w08;
static uint64_t *galois_split_w64[15] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };

uint64_t galois_shift_multiply_w64 (uint64_t x, uint64_t y)
{
  int i, k;
  uint64_t prod, j;
  uint64_t scratch[64];
  uint64_t one;
  uint64_t pp;

  one = 1;
  pp = 0;
  pp = 0x1b;

  for (i = 0; i < 64; i++) {
    scratch[i] = y;

    if (y & (one << 63)) {
      y = y << 1;
      y ^= pp;
    } else {
      y = y << 1;
    }
  }

  prod = 0;
  for (i = 0; i < 64; i++) {
    if (x & (one << i)) {
      j = 1;
      for (k = 0; k < 64; k++) {
	prod = prod ^ (j & scratch[i]);
	j <<= 1;
      }
    }
  }

  return prod;
}

uint64_t galois_shift_divide_w64 (uint64_t x, uint64_t y)
{
  uint64_t inverse;
  if (y == 0) return -1;
  if (x == 0) return 0;

  inverse = galois_shift_inverse_w64(y);
  return galois_shift_multiply_w64(x, inverse);
}

uint64_t galois_create_split_w64 ()
{
  uint64_t p1, p2, i, j, p1elt, p2elt, index, ishift, jshift, *table;
  uint64_t a, b;

  for (i = 0; i < 15; i++) {
    galois_split_w64[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 8) * (1 << 8));
    if (galois_split_w64[i] == NULL) {
      for (i--; i >= 0; i--) free(galois_split_w64[i]);
      return -1;
    }
  }

  for (i = 0; i < 8; i+=7) {
    ishift = i * 8;
    for (j = (i == 0)? 0 : 1; j < 8; j++) {
      jshift = j * 8;
      table = galois_split_w64[i+j];
      index = 0;
      for (p1 = 0; p1 < (1 << 8); p1++) {
	p1elt = (p1 << ishift);
	for (p2 = 0; p2 < (1 << 8); p2++) {
	  p2elt = (p2 << jshift);
	  table[index] = galois_shift_multiply_w64(p1elt, p2elt);
	  index++;
	}
      } 
    }
  }

}

uint64_t galois_split_multiply_w64 (uint64_t x, uint64_t y)
{
  uint64_t a, b, accumulator, i8, j8;
  uint64_t i, j;
  uint64_t index_mod;

  accumulator = 0;
  index_mod = 255;

  i8 = 0;
  for (i = 0; i < 8; i++) {
    a = (((x >> i8) & index_mod) << 8);
    j8 = 0;
    for (j = 0; j < 8; j++) {
      b = ((y >> j8) & index_mod);
      accumulator ^= galois_split_w64[i+j][a|b];
      j8 += 8;
    }
    i8 += 8;
  }
  return accumulator;

}

uint64_t galois_split_divide_w64 (uint64_t x, uint64_t y)
{
  uint64_t inverse;
  if (y == 0) return -1;
  if (x == 0) return 0;

  inverse = galois_shift_inverse_w64(y);
  return galois_shift_multiply_w64(x, inverse);
}

void galois_invert_binary_matrix_w64 (uint64_t *mat, uint64_t *inv)
{
  uint64_t cols, rows;
  int i, j, k;
  uint64_t tmp;
  uint64_t one;

  one = 1;
  cols = 64;
  rows = 64;

  for (i = 0; i < rows; i++) {
    inv[i] = 0;
    inv[i] = (one << i);
  }

  /* First -- convert into upper triangular */

  for (i = 0; i < cols; i++) {

    /* Swap rows if we ave a zero i,i element.  If we can't swap, then the 
       matrix was not invertible */

    if ((mat[i] & (one << i)) == 0) {
      for (j = i+1; 
	  j < rows && (mat[j] & (one << i)) == 0; 
	  j++);
      if (j == rows) {
	fprintf(stderr, "galois_invert_matrix: Matrix not invertible!!\n");
	exit(1);
      }
      tmp = mat[i]; mat[i] = mat[j]; mat[j] = tmp;
      tmp = inv[i]; inv[i] = inv[j]; inv[j] = tmp;
    }

    /* Now for each j>i, add A_ji*Ai to Aj */
    for (j = i+1; j != rows; j++) {
      if ((mat[j] & (one << i)) != 0) {
	mat[j] ^= mat[i];
	inv[j] ^= inv[i];
      }
    }
  }

  /* Now the matrix is upper triangular.  Start at the top and multiply down */

  for (i = rows-1; i >= 0; i--) {
    for (j = 0; j < i; j++) {
      if (mat[j] & (one << i)) {
	inv[j] ^= inv[i];
      }
    }
  }

}

uint64_t galois_shift_inverse_w64 (uint64_t x)
{
  uint64_t mat[64];
  uint64_t inv[64];
  uint64_t ind, i, j, k, prod;
  uint64_t pp, one;

  one = 1;
  pp = 0;
  pp = 0x1b;

  for (i = 0; i < 64; i++) {
    mat[i] = x;

    if (x & (one << 63)) {
      x = x << 1;
      x = x ^ pp;
    } else {
      x = x << 1;
    }
  }

  galois_invert_binary_matrix_w64(mat, inv);

  return inv[0];
}

void print_big (uint64_t a, int size)
{
  int i;
  unsigned char *b;
  b = (unsigned char *) &a;
  for (i = size-1; i >= 0; i--) {
    if (b[i/8] & (1 << i%8)) printf("1");
    else printf("0");
    if (i % 8 == 0) printf(" ");
  }
}

uint32_t galois_composite_w8_multiply_k2 (uint32_t a, uint32_t b, uint32_t n)
{
}

uint64_t galois_composite_multiply_k2 (uint64_t a, uint64_t b, uint64_t n/*, uint64_t s*/)
{
  if (a == 0 || b == 0) return 0;
  uint64_t l;
  uint64_t s;
  uint64_t product, mask0;
  uint64_t p2, p3;
	uint64_t one;

  one = 1;
  l = n/2;
	///*
  if (n == 16) s = 63;
  else if (n == 32) s = 40188;
  else if (n == 8) s = 4;
	else s = (1 << l) -1;
	//*/
  mask0 = (one << l) - 1;

  p2 = galois_single_multiply(a >> l, b >> l, l);
  p3 = p2;
  p2 = galois_single_multiply(s, p2, l);

  product = galois_single_multiply(a >> l, b & mask0, l);
  product ^= galois_single_multiply(a & mask0, b >> l, l);
  product ^= p2;
  product <<= l;
  product ^= p3;
  product ^= galois_single_multiply(a & mask0, b & mask0, l);

  return product;
}

uint64_t galois_composite_inverse_k2 (uint64_t a, uint64_t n/*, uint64_t s*/)
{
  if (a <= 0) return 0;
  uint64_t l;
  uint64_t s;
  uint64_t c;
  uint64_t mask0, inv;
  uint64_t a1i, a0i;
	uint64_t one;

	one = 1;
  l = n/2;
	///*
  if (n == 16) s = 63;
  else if (n == 32) s = 40188;
  else if (n == 8) s = 4;
	else s = (1 << l) -1;
	//*/
  mask0 = (one << l) - 1;

  if ((a >> l) & mask0 == 0) {
    inv = galois_single_divide(1, a & mask0, l);
  } else if ((a & mask0) == 0) {
    a1i = galois_single_divide(1, (a >> l) & mask0, l);
    inv = a1i;
    inv <<= l;
    inv ^= galois_single_multiply(a1i, s, l);
  } else {
    a1i = galois_single_divide(1, (a >> l) & mask0, l);
    a0i = galois_single_divide(1, a & mask0, l);

    c = galois_single_multiply((a >> l) & mask0, a0i, l);
    c ^= galois_single_multiply(a & mask0, a1i, l);
    c ^= s & mask0;
    c = galois_single_divide(1, c & mask0, l);
    c = galois_single_multiply(c & mask0, (a >> l) & mask0, l);
    c = galois_single_multiply(c & mask0, a0i, l);

    inv = galois_single_multiply(a1i, c & mask0, l);
    inv <<= l;
    inv ^= galois_single_multiply(c ^ 1, a0i, l);
  }
  return inv;
}

uint64_t galois_composite_divide_k2 (uint64_t a, uint64_t b, uint64_t n/*, uint64_t s*/)
{
  uint64_t inverse;
  if (b == 0) return -1;
  if (a == 0) return 0;

  inverse = galois_composite_inverse_k2 (b, n);
  return galois_composite_multiply_k2 (a, inverse, n);
}

void galois_composite_w08_region_mult_k2 (unsigned char *region,
    int multby,
    int nbytes,
    unsigned char *r2,
    int add)
{
  unsigned char *ur1, *ur2, *cp;
  unsigned char prod;
  uint64_t l, *lp2, *lptop;
  unsigned char *lp;
  int sol, i, j;

  int *table;
  uint32_t mask0;
  uint32_t s;
  uint32_t mb0, mb1;
  uint32_t b0, b1;
  uint32_t p2;

  ur1 = region;
  ur2 = (r2 == NULL) ? ur1 : r2;

  if (galois_create_mult_tables(4) < 0) {
    fprintf(stderr, "galois_08_region_multiply -- couldn't make multiplication tables\n");
    exit(1);
  }

  table = galois_get_mult_table(4);

  mask0 = (1 << 4) - 1;
  mb0 = (multby & mask0) << 4;
  mb1 = (multby >> 4) << 4;
  s = 4;

  if (multby == 0) {
    if (!add) {
      lp2 = (uint64_t *) ur2;
      ur2 += nbytes;
      lptop = (uint64_t *) ur2;
      while (lp2 < lptop) { *lp2 = 0; lp2++; }
    }
    return;
  }

  if (r2 == NULL || !add) {
    for (i = 0; i < nbytes; i++) {
      b0 = ur1[i] & mask0;
      b1 = ur1[i] >> 4;
      p2 = table[mb1 | b1];

      prod = table[mb1 | b0];
      prod ^= table[mb0 | b1];
      prod ^= table[(p2 << 4) | s];
      prod <<= 4;
      prod ^= p2;
      prod ^= table[mb0 | b0];

      ur2[i] = prod;
    }
    return;
  } else {
    sol = sizeof(uint64_t);
    lp2 = &l;
    lp = (unsigned char *) lp2;
    for (i = 0; i < nbytes; i += sol) {
      cp = ur2+i;
      lp2 = (uint64_t *) cp;
      for (j = 0; j < sol; j++) {
	b0 = ur1[i+j] & mask0;
	b1 = ur1[i+j] >> 4;
	p2 = table[mb1 | b1];

	prod = table[mb1 | b0];
	prod ^= table[mb0 | b1];
	prod ^= table[(p2 << 4) | s];
	prod <<= 4;
	prod ^= p2;
	prod ^= table[mb0 | b0];

	lp[j] = prod;
      }
      *lp2 = (*lp2) ^ l;
    }
  }
  return;
}

unsigned char *galois_small_table(int w) {
  int j, x, y, logx;

  if (table_w08 != NULL) return table_w08;
  table_w08 = (unsigned char *) malloc(sizeof(unsigned char) * (1 << 8) * (1 << 8));
  if (table_w08 == NULL) return NULL;

  /* Set mult/div tables for x = 0 */
  j = 0;
  table_w08[j] = 0;   /* y = 0 */
  j++;
  for (y = 1; y < (1 << 8); y++) {   /* y > 0 */
    table_w08[j] = 0;
    j++;
  }

  for (x = 1; x < (1 << 8); x++) {  /* x > 0 */
    table_w08[j] = 0; /* y = 0 */
    j++;
    for (y = 1; y < (1 << 8); y++) {  /* y > 0 */
      table_w08[j] = (unsigned char) galois_single_multiply(x, y, 8);
      j++;
    }
  }
  return table_w08;
}

void galois_composite_w16_region_mult_k2 (unsigned char *region,
    int multby,
    int nbytes,
    unsigned char *r2,
    int add)
{
  unsigned short *ur1, *ur2, *cp;
  int prod;
  uint64_t l, *lp2, *lptop;
  unsigned short *lp;
  int sol, i, j;

  unsigned char *table;
  uint32_t mask0;
  uint32_t s;
  uint32_t mb0, mb1;
  uint32_t b0, b1;
  uint32_t p2;


  ur1 = (unsigned short *) region;
  ur2 = (r2 == NULL) ? ur1 : (unsigned short *) r2;
  nbytes /= 2;

  if (multby == 0) {
    if (!add) {
      lp2 = (uint64_t *) ur2;
      ur2 += nbytes;
      lptop = (uint64_t *) ur2;
      while (lp2 < lptop) { *lp2 = 0; lp2++; }
    }
    return;
  }

  /*
     if (galois_create_mult_tables(8) < 0) {
     fprintf(stderr, "galois_16_region_multiply -- couldn't make log tables\n");
     exit(1);
     }
   */

  table = galois_small_table(8);
  mask0 = (1 << 8) - 1;
  s = 63;
  mb0 = (multby & mask0) << 8;
  mb1 = (multby >> 8) << 8;

  if (r2 == NULL || !add) {
    for (i = 0; i < nbytes; i++) {
      if (ur1[i] == 0) {
	ur2[i] = 0;
      } else {
	b0 = ur1[i] & mask0;
	b1 = ur1[i] >> 8;
	p2 = table[mb1 | b1];

	prod = table[mb1 | b0];
	prod ^= table[mb0 | b1];
	prod ^= table[(p2 << 8) | s];
	prod <<= 8;
	prod ^= p2;
	prod ^= table[mb0 | b0];

	ur2[i] = prod;
			}
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (unsigned short *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					b0 = ur1[i+j] & mask0;
					b1 = ur1[i+j] >> 8;
					p2 = table[mb1 | b1];

					prod = table[mb1 | b0];
					prod ^= table[mb0 | b1];
					prod ^= table[(p2 << 8) | s];
					prod <<= 8;
					prod ^= p2;
					prod ^= table[mb0 | b0];

					lp[j] = prod;
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

void galois_split_composite_w16_region_mult_k2 (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	unsigned char m1, m0, sm1;
	unsigned int mask0;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = &region[nbytes/2];
	dhp = (r2 == NULL) ? shp : r2;
	dlp = &dhp[nbytes/2];
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* create w08 mult tables */
	if (galois_create_mult_tables(8) < 0) {
		fprintf(stderr, "galois_16_region_multiply -- couldn't make log tables\n");
		exit(1);
	}

	/* set up constant multipliers */
	mask0 = (1 << 8) -1;
	m1 = (multby >> 8) & mask0;
	m0 = (multby & mask0);
	sm1 = galois_single_multiply(m1, 63, 8);

	if (r2 == NULL) {
		/*****  THIS DOESN'T FUNCTION FOR R2 == NULL AND ADD == 1 *****/
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = &tmp[nbytes];
		galois_w08_region_multiply((char *) slp, m1, nbytes, (char *) dhp, add);
		galois_w08_region_multiply((char *) shp, m0, nbytes, (char *) dhp, 1);
		galois_w08_region_multiply((char *) shp, sm1, nbytes, (char *) dhp, 1);

		galois_w08_region_multiply((char *) shp, m1, nbytes, (char *) dlp, add);
		galois_w08_region_multiply((char *) slp, m0, nbytes, (char *) dlp, 1);
		
		free (region);
		region = tmp;

		return;
	}

	/* w16 region multiplies */
	galois_w08_region_multiply((char *) slp, m1, nbytes, (char *) dhp, add);
	galois_w08_region_multiply((char *) shp, m0, nbytes, (char *) dhp, 1);
	galois_w08_region_multiply((char *) shp, sm1, nbytes, (char *) dhp, 1);

	galois_w08_region_multiply((char *) shp, m1, nbytes, (char *) dlp, add);
	galois_w08_region_multiply((char *) slp, m0, nbytes, (char *) dlp, 1);

	return;
}


void galois_split_composite_w32_region_mult_k2 (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* low region pointer, high region pointer for destination and source */
	/* multby split into high, low, and s * high */
	int i, j;
	unsigned char *dhp, *shp, *dlp, *slp;
	unsigned short m1, m0, sm1;
	unsigned int mask0;
	unsigned char *tmp;

	/* set region pointers */
	shp = region;
	slp = &region[nbytes/2];
	dhp = (r2 == NULL) ? shp : r2;
	dlp = &dhp[nbytes/2];
	nbytes = nbytes/2;

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dhp, nbytes);
			bzero (dlp, nbytes);
		}
		return;
	}

	/* create w16 log tables */
	if (galois_create_log_tables(16) < 0) {
		fprintf(stderr, "galois_32_region_multiply -- couldn't make log tables\n");
		exit(1);
	}

	/* set up constant multipliers */
	mask0 = (1 << 16) -1;
	m1 = (multby >> 16) & mask0;
	m0 = (multby & mask0);
	sm1 = galois_single_multiply(m1, 40188, 16);

	if (r2 == NULL) {
		tmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*2);
		dhp = tmp;
		dlp = &tmp[nbytes];
		galois_w16_region_multiply((char *) slp, m1, nbytes, (char *) dhp, add);
		galois_w16_region_multiply((char *) shp, m0, nbytes, (char *) dhp, 1);
		galois_w16_region_multiply((char *) shp, sm1, nbytes, (char *) dhp, 1);

		galois_w16_region_multiply((char *) shp, m1, nbytes, (char *) dlp, add);
		galois_w16_region_multiply((char *) slp, m0, nbytes, (char *) dlp, 1);
		
		free (region);
		region = tmp;

		return;
	}

	/* w16 region multiplies */
	galois_w16_region_multiply((char *) slp, m1, nbytes, (char *) dhp, add);
	galois_w16_region_multiply((char *) shp, m0, nbytes, (char *) dhp, 1);
	galois_w16_region_multiply((char *) shp, sm1, nbytes, (char *) dhp, 1);

	galois_w16_region_multiply((char *) shp, m1, nbytes, (char *) dlp, add);
	galois_w16_region_multiply((char *) slp, m0, nbytes, (char *) dlp, 1);

	return;
}

void galois_composite_w32_region_mult_k2 (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2, *cp, *lp, l, *lp2;
	uint32_t i, j;

	int *logtable, *ilogtable;
	uint32_t mblog1, mblog0, mb;
	uint32_t prod;
	uint32_t mask0, s, blog1, blog0, p2, p3;

	mb = multby;
	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(int);

	if ((int) multby == 0) {
		if (!add) {
			lp = (uint32_t *) ur2;
			ur2 += nbytes;
			lp2 = (uint32_t *) ur2;
			while (lp < lp2) { *lp = 0; lp++; }
		}
		return;
	}

	if (galois_create_log_tables(16) < 0) {
		fprintf(stderr, "galois_32_region_multiply -- couldn't make log tables\n");
		exit(1);
	}

	logtable = galois_get_log_table(16);
	ilogtable = galois_get_ilog_table(16);

	s = 40188;
	mask0 = (1 << 16) - 1;

	mblog1 = logtable[mb >> 16];
	mblog0 = logtable[mb & mask0];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if ((int) ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				p2 = galois_single_multiply(mb >> 16, ur1[i] >> 16, 16);
				p3 = p2;
				p2 = galois_single_multiply(p2, s, 16);

				prod = galois_single_multiply(mb >> 16, ur1[i] & mask0, 16);
				prod ^= galois_single_multiply(mb & mask0, ur1[i] >> 16, 16);
				prod ^= p2;
				prod <<= 16;
				prod ^= p3;
				prod ^= galois_single_multiply(mb & mask0, ur1[i] & mask0, 16);

				//prod = galois_composite_multiply_k2(multby, ur1[i], 32);

				/*
					 blog0 = logtable[ur1[i] & mask0];
					 blog1 = logtable[ur1[i] >> 16];
					 p2 = ilogtable[blog1 + mblog1];
					 p3 = p2;
					 p2 = ilogtable[logtable[p2] + logtable[s]];

					 prod  = ilogtable[mblog1 + blog0];
					 prod ^= ilogtable[mblog0 + blog1];
					 prod ^= p2;

					 prod <<= 16;

					 prod ^= ilogtable[mblog0 + blog0];
					 prod ^= p3;
				 */

				ur2[i] = prod;
			}
		}
	} else {
		lp = (uint32_t *) &l;
		for (i = 0; i < nbytes; i ++) {
			cp = ur2+i;
			lp2 = (uint32_t *) cp;
			if (ur1[i] == 0) {
				l = 0;
			} else {
				p2 = galois_single_multiply(mb >> 16, ur1[i] >> 16, 16);

				prod = galois_single_multiply(mb >> 16, ur1[i] & mask0, 16);
				prod ^= galois_single_multiply(mb & mask0, ur1[i] >> 16, 16);
				prod ^= galois_single_multiply(s, p2, 16);
				prod <<= 16;
				prod ^= p2;
				prod ^= galois_single_multiply(mb & mask0, ur1[i] & mask0, 16);

				/*
					 blog0 = logtable[ur1[i] & mask0];
					 blog1 = logtable[ur1[i] >> 16];
					 p2 = ilogtable[blog1 + mblog1];

					 prod  = ilogtable[mblog1 + blog0];
					 prod ^= ilogtable[mblog0 + blog1];
					 prod ^= ilogtable[logtable[p2] + logtable[s]];

					 prod <<= 16;

					 prod ^= p2;
					 prod ^= ilogtable[mblog0 + blog0];
				 */
				//prod = galois_composite_multiply_k2(multby, ur1[i], 32);

				l = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}


uint64_t galois_composite_multiply_k4 (uint64_t a, uint64_t b, uint64_t n)
{
	uint64_t l;
	uint64_t s, t;
	uint64_t product, mask1, mask0, mask2, mask3;
	uint64_t a4, a5, a6;
	uint64_t sa4, ta4, sa5, ta5, sa6, ta6;

	l = n/4;
	s = 3;
	t = 8;

	mask0 = (1 << l) - 1;
	mask1 = mask0 << l;
	mask2 = mask1 << l;
	mask3 = mask2 << l;

	a4  = galois_single_multiply((a & mask1) >> (l * 1), (b & mask3) >> (l * 3), l);
	a4 ^= galois_single_multiply((a & mask2) >> (l * 2), (b & mask2) >> (l * 2), l);
	a4 ^= galois_single_multiply((a & mask3) >> (l * 3), (b & mask1) >> (l * 1), l);

	a5  = galois_single_multiply((a & mask3) >> (l * 3), (b & mask2) >> (l * 2), l);
	a5 ^= galois_single_multiply((a & mask2) >> (l * 2), (b & mask3) >> (l * 3), l);

	a6  = galois_single_multiply((a & mask3) >> (l * 3), (b & mask3) >> (l * 3), l);

	sa4 = galois_single_multiply(s, a4, l);
	ta4 = galois_single_multiply(t, a4, l);
	sa5 = galois_single_multiply(s, a5, l);
	ta5 = galois_single_multiply(t, a5, l);
	sa6 = galois_single_multiply(s, a6, l);
	ta6 = galois_single_multiply(t, a6, l);

	/* x ^ 3 terms */
	product  = galois_single_multiply( a          >> (l * 3), (b & mask0),            l);
	product ^= galois_single_multiply((a & mask2) >> (l * 2), (b & mask1) >> (l * 1), l);
	product ^= galois_single_multiply((a & mask1) >> (l * 1), (b & mask2) >> (l * 2), l);
	product ^= galois_single_multiply((a & mask0),             b          >> (l * 3), l);

	product ^= a5;
	product ^= sa6;
	product <<= l;

	/* x ^ 2 terms */
	product ^= galois_single_multiply((a & mask2) >> (l * 2), (b & mask0), l);
	product ^= galois_single_multiply((a & mask1) >> (l * 1), (b & mask1) >> (l * 1), l);
	product ^= galois_single_multiply((a & mask0),            (b & mask2) >> (l * 2), l);

	product ^= a4;
	product ^= sa5;
	product ^= ta6;
	product ^= a6;
	product <<= l;

	/* x ^ 1 terms */
	product ^= galois_single_multiply((a & mask1) >> (l * 1), (b & mask0), l);
	product ^= galois_single_multiply((a & mask0),            (b & mask1) >> (l * 1), l);
	product ^= sa4;
	product ^= ta5;
	product ^= sa6;
	product <<= l;

	/* x ^ 0 terms */
	product ^= galois_single_multiply((a & mask0),            (b & mask0), l);
	product ^= ta4;
	product ^= ta6;

	return product;
}

uint64_t galois_composite_divide_k4 (uint64_t a, uint64_t b, uint64_t n)
{
	uint64_t inverse;
	if (b == 0) return -1;
	if (a == 0) return 0;
	inverse = galois_composite_inverse_k4(b, n);
	return galois_composite_multiply_k4(inverse, a, n);
}

uint64_t galois_composite_inverse_k4 (uint64_t a, uint64_t n)
{
	uint32_t remainder[2][5], r1[5];
	uint32_t a1[5], a0[5], tmp[6];
	uint32_t quotient[5];
	uint32_t inv;

	int mask0;
	int s, t;
	int high_bit, high_bit_r, shift_count, i, j;

	mask0 = (1 << (n/4)) -1;
	s = 3;
	t = 8;

	remainder[0][0] = t;
	remainder[0][1] = s;
	remainder[0][2] = 1;
	remainder[0][3] = 0;
	remainder[0][4] = 1;

	remainder[1][0] = a & mask0;
	remainder[1][1] = (a >> (n/4)) & mask0;
	remainder[1][2] = (a >> (2*(n/4))) & mask0;
	remainder[1][3] = (a >> (3*(n/4))) & mask0;
	remainder[1][4] = 0;

	if (a == 0) return -1;
	for (high_bit_r = 4; remainder[1][high_bit_r] == 0; high_bit_r--);
	if (high_bit_r == 0) return galois_inverse(remainder[1][0], (n/4));

	for (i = 0; i < 5; i++) a1[i] = a0[i] = 0;
	a1[0] = 1;
	high_bit = 4;
	high_bit_r = 4;
	while (high_bit > 0) {
		for (i = 0; i < 5; i++) quotient[i] = 0;
		high_bit = high_bit_r;

		for (high_bit_r = 4; remainder[1][high_bit_r] == 0; high_bit_r--);

		/* This is the polynomial division part */
		while(high_bit >= high_bit_r) {
			for (i = 0; i < 5; i++) r1[i] = remainder[1][i];
			shift_count = 0;

			shift_count = high_bit - high_bit_r;
			inv = galois_single_divide(remainder[0][high_bit], r1[high_bit_r], (n/4));
			for (i = 0; i < 5; i++) r1[i] = galois_single_multiply(inv, r1[i], (n/4));

			/* set the quotient */
			quotient[shift_count] = inv;

			/* next remainder in the division process */
			for (i = 0; i <= high_bit_r; i++) remainder[0][i+shift_count] ^= r1[i];
			while(remainder[0][high_bit] == 0) high_bit--;
		}

		for (i = 4; i >= 0; i--) r1[i] = remainder[0][i];
		for (i = 4; i >= 0; i--) remainder[0][i] = remainder[1][i];
		for (i = 4; i >= 0; i--) remainder[1][i] = r1[i]; 
		for (i = 5; i >= 0; i--) tmp[i] = 0;
		for (i = 3; i >= 0; i--) {
			if (a1[i] > 0) {
				for (j = 3; j >= 0; j--) {
					if (quotient[j] > 0) {
						tmp[i+j] = tmp[i+j] ^ (galois_single_multiply(quotient[j], a1[i], (n/4)));
					}
				}
			}
		}
		for (i = 3; i >= 0; i--) { tmp[i] ^= a0[i]; a0[i] = a1[i]; }
		for (i = 3; i >= 0; i--) a1[i] = tmp[i];
	}
	if (remainder[1][0] > 1) {
		inv = galois_inverse(remainder[1][0], (n/4));
		for (i = 0; i < 4; i++) {
			a1[i] = galois_single_multiply(inv, a1[i], (n/4));
		}
	}
	inv = 0;
	for (i = 3; i >= 0; i--) { inv <<= (n/4); inv ^= a1[i]; }
	return inv;
}

void galois_split_composite_w32_region_mult_k4 (unsigned char *region,
		uint32_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	/* destination and source pointers: 3 is most significant, 0 is least */
	/* multby split into four, and s && t */
	int i, j;
	unsigned char *dp3, *dp2, *dp1, *dp0, *sp3, *sp2, *sp1, *sp0;
	uint32_t m3, m2, m1, m0;
	uint32_t sm3, sm2, sm1, s;
	uint32_t tm3, tm2, tm1, tp1m3, t;
	uint32_t mask0;
	unsigned char *rtmp;

	/* set region pointers */
	nbytes = nbytes/4;
	sp3 = region;
	sp2 = &sp3[nbytes];
	sp1 = &sp2[nbytes];
	sp0 = &sp1[nbytes];

	dp3 = (r2 == NULL) ? sp3: r2;
	dp2 = &dp3[nbytes];
	dp1 = &dp2[nbytes];
	dp0 = &dp1[nbytes];

	/* is multby 0? Speed that up */
	if (multby == 0) {
		if (!add) {
			bzero (dp3, nbytes);
			bzero (dp2, nbytes);
			bzero (dp1, nbytes);
			bzero (dp0, nbytes);
		}
		return;
	}

	/* create w08 mult tables */
	if (galois_create_mult_tables(8) < 0) {
		fprintf(stderr, "galois_32_region_multiply -- couldn't make mult tables\n");
		exit(1);
	}

	/* set up constant multipliers */
	s = 3;
	t = 8;
	mask0 = (1 << 8) -1;
	m3 = (multby >> 24) & mask0;
	m2 = (multby >> 16) & mask0;
	m1 = (multby >> 8) & mask0;
	m0 = (multby & mask0);

	sm3 = galois_single_multiply((int) s, (int) m3, 8);
	sm2 = galois_single_multiply((int) s, (int) m2, 8);
	sm1 = galois_single_multiply((int) s, (int) m1, 8);

	tm3 = galois_single_multiply((int) t, (int) m3, 8);
	tm2 = galois_single_multiply((int) t, (int) m2, 8);
	tm1 = galois_single_multiply((int) t, (int) m1, 8);
	tp1m3 = galois_single_multiply((int) t ^ 1, (int) m3, 8);

	if (r2 == NULL) {
		rtmp = (unsigned char *) malloc(sizeof(unsigned char) * nbytes*4);
		dp3 = rtmp;
		dp2 = &dp3[nbytes];
		dp1 = &dp2[nbytes];
		dp0 = &dp1[nbytes];
	}

	galois_w08_region_multiply((char *) sp0, m3, nbytes, (char *) dp3, add);
	galois_w08_region_multiply((char *) sp1, m2, nbytes, (char *) dp3, 1);
	galois_w08_region_multiply((char *) sp2, m1, nbytes, (char *) dp3, 1);
	galois_w08_region_multiply((char *) sp2, m3, nbytes, (char *) dp3, 1);
	galois_w08_region_multiply((char *) sp3, m0, nbytes, (char *) dp3, 1);
	galois_w08_region_multiply((char *) sp3, m2, nbytes, (char *) dp3, 1);
	galois_w08_region_multiply((char *) sp3, sm3, nbytes, (char *) dp3, 1);

	galois_w08_region_multiply((char *) sp0, m2, nbytes, (char *) dp2, add);
	galois_w08_region_multiply((char *) sp1, m1, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp1, m3, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp2, m0, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp2, m2, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp2, sm3, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp3, m1, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp3, sm2, nbytes, (char *) dp2, 1);
	galois_w08_region_multiply((char *) sp3, tp1m3, nbytes, (char *) dp2, 1);

	galois_w08_region_multiply((char *) sp0, m1, nbytes, (char *) dp1, add);
	galois_w08_region_multiply((char *) sp1, m0, nbytes, (char *) dp1, 1);
	galois_w08_region_multiply((char *) sp1, sm3, nbytes, (char *) dp1, 1);
	galois_w08_region_multiply((char *) sp2, sm2, nbytes, (char *) dp1, 1);
	galois_w08_region_multiply((char *) sp2, tm3, nbytes, (char *) dp1, 1);
	galois_w08_region_multiply((char *) sp3, sm1, nbytes, (char *) dp1, 1);
	galois_w08_region_multiply((char *) sp3, tm2, nbytes, (char *) dp1, 1);
	galois_w08_region_multiply((char *) sp3, sm3, nbytes, (char *) dp1, 1);

	galois_w08_region_multiply((char *) sp0, m0, nbytes, (char *) dp0, add);
	galois_w08_region_multiply((char *) sp1, tm3, nbytes, (char *) dp0, 1);
	galois_w08_region_multiply((char *) sp2, tm2, nbytes, (char *) dp0, 1);
	galois_w08_region_multiply((char *) sp3, tm1, nbytes, (char *) dp0, 1);
	galois_w08_region_multiply((char *) sp3, tm3, nbytes, (char *) dp0, 1);

	if (r2 == NULL) { free(region); region = rtmp; }
	return;

}

void galois_composite_w32_region_mult_k4 (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2, *cp, *lp, l, *lp2;
	uint32_t prod;
	int sol, i, j;

	unsigned char *table;
	uint32_t mask0, mask1, mask2, mask3;
	uint32_t s, t;
	uint32_t mb0, mb1, mb2, mb3;
	uint32_t b0, b1, b2, b3;
	uint32_t a4, a5, a6;
	uint32_t sa4, ta4, sa5, ta5, sa6, ta6;

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);

	if (multby == 0) {
		if (!add) {
			lp = (uint32_t *) ur2;
			ur2 += nbytes;
			lp2 = (uint32_t *) ur2;
			while(lp < lp2) { *lp = 0; lp++; }
		}
		return;
	}

	table = galois_small_table(8);
	s = 3;
	t = 8;

	mask0 = (1 << 8) - 1;
	mask1 = mask0 << 8;
	mask2 = mask1 << 8;
	mask3 = mask2 << 8;

	mb0 = (multby & mask0) << 8;
	mb1 = ((multby & mask1) >> 8) << 8;
	mb2 = ((multby & mask2) >> 16) << 8;
	mb3 = ((multby & mask3) >> 24) << 8;

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				b0 = ur1[i] & mask0;
				b1 = (ur1[i] & mask1) >> 8;
				b2 = (ur1[i] & mask2) >> 16;
				b3 = (ur1[i] & mask3) >> 24;

				a4  = table[mb1 | b3];
				a4 ^= table[mb2 | b2];
				if (mb3 < 0) printf("%d\n", mb3);
				a4 ^= table[mb3 | b1];

				a5  = table[mb3 | b2];
				a5 ^= table[mb2 | b3];

				a6  = table[mb3 | b3];

				sa4 = table[(a4 << 8) | s];
				ta4 = table[(a4 << 8) | t];
				sa5 = table[(a5 << 8) | s];
				ta5 = table[(a5 << 8) | t];
				sa6 = table[(a6 << 8) | s];
				ta6 = table[(a6 << 8) | t];

				/* x ^ 3 terms */
				prod  = table[mb3 | b0];
				prod ^= table[mb2 | b1];
				prod ^= table[mb1 | b2];
				prod ^= table[mb0 | b3];
				prod ^= a5;
				prod ^= sa6;
				prod <<= 8;

				/* x ^ 2 terms */
				prod ^= table[mb2 | b0];
				prod ^= table[mb1 | b1];
				prod ^= table[mb0 | b2];
				prod ^= a4;
				prod ^= sa5;
				prod ^= ta6;
				prod ^= a6;
				prod <<= 8;

				/* x ^ 1 terms */
				prod ^= table[mb1 | b0];
				prod ^= table[mb0 | b1];
				prod ^= sa4;
				prod ^= ta5;
				prod ^= sa6;
				prod <<= 8;

				/* x ^ 0 terms */
				prod ^= table[mb0 | b0];
				prod ^= ta4;
				prod ^= ta6;

				ur2[i] = prod;
			}
		}
	} else {
		lp = (uint32_t *) &l;
		for (i = 0; i < nbytes; i++) {
			cp = ur2 + i;
			lp2 = (uint32_t *) cp;
			if (ur1[i] == 0) {
				l = 0;
			} else {
				b0 = ur1[i] & mask0;
				b1 = (ur1[i] & mask1) >> 8;
				b2 = (ur1[i] & mask2) >> 16;
				b3 = (ur1[i] & mask3) >> 24;

				a4  = table[mb1 | b3];
				a4 ^= table[mb2 | b2];
				a4 ^= table[mb3 | b1];

				a5  = table[mb3 | b2];
				a5 ^= table[mb2 | b3];

				a6  = table[mb3 | b3];

				sa4 = table[(a4 << 8) | s];
				ta4 = table[(a4 << 8) | t];
				sa5 = table[(a5 << 8) | s];
				ta5 = table[(a5 << 8) | t];
				sa6 = table[(a6 << 8) | s];
				ta6 = table[(a6 << 8) | t];

				/* x ^ 3 terms */
				prod  = table[mb3 | b0];
				prod ^= table[mb2 | b1];
				prod ^= table[mb1 | b2];
				prod ^= table[mb0 | b3];
				prod ^= a5;
				prod ^= sa6;
				prod <<= 8;

				/* x ^ 2 terms */
				prod ^= table[mb2 | b0];
				prod ^= table[mb1 | b1];
				prod ^= table[mb0 | b2];
				prod ^= a4;
				prod ^= sa5;
				prod ^= ta6;
				prod ^= a6;
				prod <<= 8;

				/* x ^ 1 terms */
				prod ^= table[mb1 | b0];
				prod ^= table[mb0 | b1];
				prod ^= sa4;
				prod ^= ta5;
				prod ^= sa6;
				prod <<= 8;

				/* x ^ 0 terms */
				prod ^= table[mb0 | b0];
				prod ^= ta4;
				prod ^= ta6;

				l = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}

	return;
}

uint64_t *galois_shift_multiply_w128 (uint64_t *x, uint64_t *y)
{
	int i, k;
	uint64_t *prod;
	uint64_t scratch[128][2];
	uint64_t *s, one, pp, j;
	uint64_t savey[2];

	prod = (uint64_t *) malloc (sizeof(uint64_t) * 2);
	one = 1;
	pp = 0x87;
	savey[0] = y[0]; savey[1] = y[1];

	for (i = 0; i < 128; i++) {
		memcpy(scratch[i], savey, sizeof(uint64_t) * 2);

		if (savey[0] & (one << 63)) {
			savey[0] <<= 1;
			if (savey[1] & (one << 63)) savey[0] ^= 1;
			savey[1] <<= 1;
			savey[1] ^= pp;
		} else {
			savey[0] <<= 1;
			if (savey[1] & (one << 63)) savey[0] ^= 1;
			y[1] <<= 1;
		}
	}

	prod[0] = 0; prod[1] = 0;
	/* Iterate over the first 64 bits of x */
	for (i = 0; i < 64; i++) {
		if (y[1] & (one << i)) {
			prod[1] = prod[1] ^ scratch[i][1];
			prod[0] = prod[0] ^ scratch[i][0];
		}
	}

	/* Iterate over the last 64 bits of x*/
	for (i = 0; i < 64; i++) {
		if (y[0] & (one << i)) {
			prod[1] = prod[1] ^ scratch[i+64][1];
			prod[0] = prod[0] ^ scratch[i+64][0];
		}
	}
	return prod;
}

uint64_t *galois_shift_divide_w128 (uint64_t *x, uint64_t *y)
{
	uint64_t *inverse;
	uint64_t *retval;
	retval = malloc(sizeof(uint64_t) * 2);
	retval[0] = 0;
	if (y[0] == 0 && y[1] == 0) { retval[1] = -1; return retval; }
	if (x[0] == 0 && x[1] == 0) { retval[1] = 0;  return retval; }
	free(retval);

	inverse = galois_shift_inverse_w128(y);
	return galois_shift_multiply_w128(x, y);
}

uint64_t *galois_shift_inverse_w128 (uint64_t *x)
{
	uint64_t remainder[2][2], auxiliary[2][2];
	uint64_t r1[2], q1[2], a1[2];
	uint64_t *tmp;
	uint64_t one;
	uint64_t *retval;
	int high_bit, high_bit_r, shift_count, i, flag;

	one = 1;
	remainder[0][0] = 0; remainder[0][1] = 0x87;

	remainder[1][0] = x[0];
	remainder[1][1] = x[1];

	auxiliary[0][0] = 0; auxiliary[0][1] = 0;
	auxiliary[0][0] = 0; auxiliary[0][1] = 1;

	high_bit   = 128;
	high_bit_r = 128;
	flag = 1;

	while (1) {
		if (remainder[1][0] == 0 && remainder[1][1] == 1) {
			break;
		}

		q1[0] = 0; q1[1] = 0;
		high_bit = high_bit_r;
		while (high_bit_r >= 64 && (remainder[1][0] & (one << (high_bit_r - 64))) == 0) high_bit_r--;
		while (high_bit_r < 64  && (remainder[1][1] & (one << high_bit_r)) == 0) high_bit_r--;

		while (high_bit >= high_bit_r) {
			r1[0] = remainder[1][0]; r1[1] = remainder[1][1];

			shift_count = high_bit - high_bit_r;
			r1[0] <<= shift_count;
			r1[0] = r1[0] ^ (r1[1] >> (64-shift_count));
			r1[1] <<= shift_count;

			if (shift_count < 64) q1[1] ^= (one << shift_count);
			else q1[0] ^= (one << (shift_count-64));

			remainder[0][0] ^= r1[0];
			remainder[0][1] ^= r1[1];
			while (high_bit >= 64 && (remainder[0][0] & (one << (high_bit - 64))) == 0) high_bit--;
			while (high_bit < 64  && (remainder[0][1] & (one << high_bit)) == 0) high_bit--;
		}

		r1[0] = remainder[0][0]; r1[1] = remainder[0][1];
		remainder[0][0] = remainder[1][0]; remainder[0][1] = remainder[1][1];
		remainder[1][0] = r1[0]; remainder[1][1] = r1[1];
		retval = galois_shift_multiply_w128(q1, auxiliary[1]);
		a1[0] = retval[0];
		a1[1] = retval[1];

		a1[0] ^= auxiliary[0][0];
		a1[1] ^= auxiliary[0][1];
		auxiliary[0][0] = auxiliary[1][0]; auxiliary[0][1] = auxiliary[1][1];
		auxiliary[1][0] = a1[0]; auxiliary[1][1] = a1[1];
	}
	retval = (uint64_t *) malloc (sizeof(uint64_t) * 2);
	bzero(retval, 16);
	memcpy(retval, &auxiliary[1], 16);
	return retval;
}

int galois_create_split_w128 ()
{
}

unsigned char *galois_split_multiply_w128 (unsigned char *x, unsigned char *y)
{
}

unsigned char *galois_split_divide_w128 (unsigned char *x, unsigned char *y)
{
}

void galois_w08_log_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *ur1, *ur2, *cp;
	int prod;
	int i, log1, j, log2;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;
	int *logtable, *ilogtable;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	if (galois_create_log_tables(8) < 0) {
		fprintf(stderr, "galois_08_region_multiply -- couldn't make log tables\n");
		exit(1);
	}
	logtable = galois_get_log_table(8);
	ilogtable = galois_get_ilog_table(8);
	log1 = logtable[multby];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				prod = logtable[ur1[i]] + log1;
				ur2[i] = ilogtable[prod];
			}
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					log2 = logtable[ur1[i+j]];
					prod = log2 + log1;
					lp[j] = ilogtable[prod];
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
}

unsigned char *galois_create_w08_lazytable_row (int row_index)
{
	unsigned char *row_elts;
	int i;
	int logx;
	int *log_table, *ilog_table;

	if (galois_create_log_tables(8) < 0) { return NULL; }

	log_table = galois_get_log_table(8);
	ilog_table = galois_get_ilog_table(8);

	row_elts = (unsigned char *) malloc(sizeof(unsigned char) * (1 << 8));

	logx = log_table[row_index];
  row_elts[0] = 0;
	for (i = 1; i < (1 << 8); i++) {
		row_elts[i] = (unsigned char) ilog_table[logx + log_table[i]];
	}
	return row_elts;
}

void galois_w08_lazytable_region_multiply (unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned char *table_row;
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		bzero(ur2, nbytes);
		return;
	}

	table_row = galois_create_w08_lazytable_row(multby);
	if (table_row == NULL) {
		fprintf(stderr, "galois_08_lazytable_region_multiply -- couldn't make multiplication tables\n");
		exit(1);
	}

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = table_row[ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = table_row[ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	free (table_row);
	return;
}

unsigned short *galois_create_w16_lazytable_row (int row_index)
{
	unsigned short *row_elts;
	int i;
	int logx;
	int *log_table, *ilog_table;

	if (galois_create_log_tables(16) < 0) { return NULL; }

	log_table = galois_get_log_table(16);
	ilog_table = galois_get_ilog_table(16);

	row_elts = (unsigned short *) malloc(sizeof(unsigned short) * (1 << 16));

	logx = log_table[row_index];
	row_elts[0] = 0;
	for (i = 1; i < (1 << 16); i++) {
		row_elts[i] = (unsigned short) ilog_table[logx + log_table[i]];
	}
	return row_elts;
}

void galois_w16_lazytable_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	unsigned short *table_row;
	unsigned short *ur1, *ur2, *cp;
	unsigned short prod;
	int i, srow, j;
	uint64_t l, *lp2;
	unsigned short *lp;
	int sol;

	ur1 = (unsigned short *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned short *) r2;
	nbytes /= 2;

	if (multby == 0) {
		bzero(ur2, nbytes*2);
		return;
	}

	table_row = galois_create_w16_lazytable_row(multby);
	if (table_row == NULL) {
		fprintf(stderr, "galois_16_lazytable_region_multiply -- couldn't make multiplication tables\n");
		exit(1);
	}

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			ur2[i] = table_row[ur1[i]];
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (unsigned short *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				lp[j] = table_row[ur1[i+j]];
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	free (table_row);
	return;
}

uint32_t **galois_w32_create_lazy_split_row(int multby)
{
	int p1, p2, index, ishift, jshift, i, j, i8, j8, k;
	uint32_t p1elt, p2elt, acache[4], **table_row;
	
	table_row = (uint32_t **) malloc(sizeof(uint32_t *) * 16);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 16; i++) {
		table_row[i] = (uint32_t *) malloc(sizeof(uint32_t) * (1 << 8));
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			return NULL;
		}
	}

	i8 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (multby >> i8) & 255;
		i8 += 8;
	}

	for (i = 0; i < 4; i++) {
		if (i == 3) ishift = i * 8;
		else ishift = 0;
		for (j = 0; j < 4; j++) {
			if (i+j > 3) {
				jshift = (i+j-3) * 8;
				ishift = 3 * 8;
			} else {
				jshift = (i+j) * 8;
				ishift = 0;
			}

			index = 0;
			p1elt = (acache[i] << ishift);
			for (p2 = 0; p2 < 256; p2++) {
				p2elt = (p2 << jshift);
				table_row[i*4+j][index] = (uint32_t) galois_shift_multiply(p1elt, p2elt, 32);
				index++;
			}
		}
	}
	return table_row;
}

void galois_w32_lazy_split_region_multiply(unsigned char *region,
		int multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2;
	int i, j, a, b, accumulator, i8, j8, k;
	int p1, p2, index, ishift, jshift;
	uint32_t p1elt, p2elt, acache[4], *table;

	uint32_t **table_row;

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint32_t));
		}
		return;
	}

/*
	i8 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (multby >> i8) & 255;
		i8 += 8;
	}
	*/

	table_row = galois_w32_create_lazy_split_row(multby);

	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*4 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		/* this implementation doesn't perform as well, but I don't know why
		for (i = 0; i < 4; i++) {
			j8 = 0;
			for (j = 0; j < 4; j++) {
				for (k = 0; k < nbytes; k++) {
					b = ((ur1[k] >> j8) & 255);
					ur2[k] = ur2[k] ^ table_row[i*4 + j][b];
				}
				j8 += 8;
			}
		}
		*/
		/*
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*4 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
		*/
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			j8 = 0;
			for (j = 0; j < 4; j++) {
				b = (ur1[k] >> j8) & 255;
				for (i = 0; i < 4; i++) {
					accumulator ^= table_row[i*4 + j][b];
				}
				j8 += 8;
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	for (i = 0; i < 16; i++) free(table_row[i]);
	free(table_row);
	return;
}

uint64_t **galois_w64_create_lazy_split_row(uint64_t multby)
{
	uint64_t p1, p2, index, ishift, jshift, i, j, i8, j8, k;
	uint64_t p1elt, p2elt, acache[8], **table_row;

	table_row = (uint64_t **) malloc(sizeof(uint64_t *) * 64);
	if (table_row == NULL) { return NULL; }

	for (i = 0; i < 64; i++) {
		table_row[i] = (uint64_t *) malloc(sizeof(uint64_t) * (1 << 8));
		if (table_row[i] == NULL) {
			for (i--; i >= 0; i--) free(table_row[i]);
			return NULL;
		}
	}

	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = (multby >> i8) & 255;
		i8 += 8;
	}

	for (i = 0; i < 8; i++) {
		for (j = 0; j < 8; j++) {
			if (i+j > 7) {
				jshift = (i+j-7) * 8;
				ishift = 7 * 8;
			} else {
				jshift = (i+j) * 8;
				ishift = 0;
			}

			index = 0;
			p1elt = (acache[i] << ishift);
			for (p2 = 0; p2 < 256; p2++) {
				p2elt = (p2 << jshift);
				table_row[i*8+j][index] = (uint64_t) galois_shift_multiply_w64(p1elt, p2elt);
				index++;
			}
		}
	}
	return table_row;
}


void galois_w64_lazy_split_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *cp, *ur2top;
	uint64_t *lp2, *lptop;
	uint64_t i, j, a, b, accumulator, i8, j8, k;
	uint64_t acache[8];
	uint64_t **table_row;

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (multby == 0) {
		if (!add) {
			bzero(r2, nbytes * sizeof(uint64_t));
		}
		return;
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	table_row = galois_w64_create_lazy_split_row(multby);
	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = (((multby >> i8) & 255) << 8);
		i8 += 8;
	}
	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*8 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= table_row[i*8 + j][b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	for (i = 0; i < 16; i++) free(table_row[i]);
	free(table_row);
	return;
}

void galois_w64_split_region_multiply (unsigned char *region,
		uint64_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint64_t *ur1, *ur2, *cp, *ur2top;
	uint64_t *lp2, *lptop;
	uint64_t i, j, a, b, accumulator, i8, j8, k;
	uint64_t acache[8];

	ur1 = (uint64_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint64_t *) r2;
	nbytes /= sizeof(uint64_t);
	ur2top = ur2 + nbytes;

	if (galois_split_w64[0]== NULL) {
		if (galois_create_split_w64() < 0) {
			fprintf(stderr, "galois_64_region_multiply -- couldn't make split multiplication tables\n"
					);
			exit(1);
		}
	}

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
		 We'll inline it here to save on the procedure call overhead */

	i8 = 0;
	for (i = 0; i < 8; i++) {
		acache[i] = (((multby >> i8) & 255) << 8);
		i8 += 8;
	}
	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_split_w64[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 8; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 8; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_split_w64[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	return;
}

void galois_w32_composite_region_mult_k22 (unsigned char *region,
		uint32_t multby,
		int nbytes,
		unsigned char *r2,
		int add)
{
	uint32_t *ur1, *ur2, *cp, *lp, l, *lp2;
	uint32_t i, j;

	uint32_t mb1, mb0, mb;
	uint32_t prod;
	uint32_t mask0, s, blog1, blog0, p2, p3;

	mb = multby;
	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);

	if ((int) multby == 0) {
		if (!add) {
			lp = (uint32_t *) ur2;
			ur2 += nbytes;
			lp2 = (uint32_t *) ur2;
			while (lp < lp2) { *lp = 0; lp++; }
		}
		return;
	}

	if (galois_create_log_tables(8) < 0) {
		fprintf(stderr, "galois_32_region_multiply -- couldn't make log tables\n");
		exit(1);
	}

	s = 40188;
	mask0 = (1 << 16) - 1;
	mb1 = (mb >> 16) & mask0;
	mb0 = (mb & mask0);

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if ((int) ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				p2 = galois_composite_multiply_k2(mb1, ur1[i] >> 16, 16);
				p3 = p2;
				p2 = galois_composite_multiply_k2(p2, s, 16);

				prod = galois_composite_multiply_k2(mb >> 16, ur1[i] & mask0, 16);
				prod ^= galois_composite_multiply_k2(mb & mask0, ur1[i] >> 16, 16);
				prod ^= p2;
				prod <<= 16;
				prod ^= p3;
				prod ^= galois_composite_multiply_k2(mb & mask0, ur1[i] & mask0, 16);

				ur2[i] = prod;
			}
		}
	} else {
		lp = (uint32_t *) &l;
		for (i = 0; i < nbytes; i ++) {
			cp = ur2+i;
			lp2 = (uint32_t *) cp;
			if (ur1[i] == 0) {
				l = 0;
			} else {
				p2 = galois_composite_multiply_k2(mb >> 16, ur1[i] >> 16, 16);

				prod = galois_composite_multiply_k2(mb >> 16, ur1[i] & mask0, 16);
				prod ^= galois_composite_multiply_k2(mb & mask0, ur1[i] >> 16, 16);
				prod ^= galois_composite_multiply_k2(s, p2, 16);
				prod <<= 16;
				prod ^= p2;
				prod ^= galois_composite_multiply_k2(mb & mask0, ur1[i] & mask0, 16);

				l = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;

}

uint64_t galois_composite_multiply_k22 (uint64_t a, uint64_t b, uint64_t n)
{
	if (a == 0 || b == 0) return 0;
	uint64_t l;
	uint64_t s;
	uint64_t product, mask0;
	uint64_t p2, p3;
	uint64_t one;

	one = 1;
	l = n/2;
	s = 40188;

	mask0 = (one << l) - 1;

	p2 = galois_composite_multiply_k2(a >> l, b >> l, l);
	p3 = p2;
	p2 = galois_composite_multiply_k2(s, p2, l);

	product = galois_composite_multiply_k2(a >> l, b & mask0, l);
	product ^= galois_composite_multiply_k2(a & mask0, b >> l, l);
	product ^= p2;
	product <<= l;
	product ^= p3;
	product ^= galois_composite_multiply_k2(a & mask0, b & mask0, l);

	return product;

}
uint64_t galois_composite_divide_k22 (uint64_t a, uint64_t b, uint64_t n)
{
	uint64_t inverse;
	if (b == 0) return -1;
	if (a == 0) return 0;

	inverse = galois_composite_inverse_k22 (b, n);
	return galois_composite_multiply_k22 (a, inverse, n);

}
uint64_t galois_composite_inverse_k22 (uint64_t a, uint64_t n)
{
	if (a <= 0) return 0;
	uint64_t l;
	uint64_t s;
	uint64_t c;
	uint64_t mask0, inv;
	uint64_t a1i, a0i;
	uint64_t one;

	one = 1;
	l = n/2;
	s = 40188;
	mask0 = (one << l) - 1;

	if ((a >> l) & mask0 == 0) {
		inv = galois_composite_divide_k2(1, a & mask0, l);
	} else if ((a & mask0) == 0) {
		a1i = galois_composite_divide_k2(1, (a >> l) & mask0, l);
		inv = a1i;
		inv <<= l;
		inv ^= galois_composite_multiply_k2(a1i, s, l);
	} else {
		a1i = galois_composite_divide_k2(1, (a >> l) & mask0, l);
		a0i = galois_composite_divide_k2(1, a & mask0, l);

		c = galois_composite_multiply_k2((a >> l) & mask0, a0i, l);
		c ^= galois_composite_multiply_k2(a & mask0, a1i, l);
		c ^= s & mask0;
		c = galois_composite_divide_k2(1, c & mask0, l);
		c = galois_composite_multiply_k2(c & mask0, (a >> l) & mask0, l);
		c = galois_composite_multiply_k2(c & mask0, a0i, l);

		inv = galois_composite_multiply_k2(a1i, c & mask0, l);
		inv <<= l;
		inv ^= galois_composite_multiply_k2(c ^ 1, a0i, l);
	}
	return inv;

}

